function [feaScore, feaIdx] = fs_unsup_OurMethod(X, k)

[nSmp, nDim] = size(X);

Xnorm = sum(X.^2, 2);
X = bsxfun(@rdivide, X, sqrt(Xnorm));

G = computeGraphViaLocalRegression(X, k);

Gs = cell(nDim, 1);
H = zeros(nDim);
f = zeros(nDim, 1);
for i1 = 1:nDim;
    Gi = computeGraphViaLocalRegression(X(:,i1), k);
    f(i1) = sum(sum(G .* Gi));
	Gs{i1} = sparse(Gi);
end

for i1 = 1:nDim
    Gi = Gs{i1};
    for i2 = 1:nDim
	    Gi2 = Gs{i2};
	    H(i1, i2) = sum(sum(Gi .* Gi2));
	end
end

H = max(H, H');
opts = [];
% opts.Display = 'iter';
opts.Display = 'off';
[feaScore, fval] = quadprog(H, -f, [], [], ones(1, nDim), 1, zeros(nDim,1), ones(nDim,1), [], opts);

[feaScore, feaIdx] = sort(feaScore, 'descend');
end

function G = computeGraphViaLocalRegression(X, k)
[nSmp, nDim] = size(X);
%****************************************************
% Compute the affinity graph via local regression
%****************************************************
Xnorm = sum(X.^2, 2);
KX = X * X'; % linear kernel
KX = (KX + KX')/2;
XD = bsxfun(@plus, Xnorm, Xnorm') - 2 * KX;
[~, Idx] = sort(XD, 2, 'ascend');
Idx = Idx(:, 2:k+1);
A = zeros(nSmp);
A(sub2ind([nSmp, nSmp], repmat([1:nSmp]', k, 1), Idx(:))) = 1;
A = max(A, A');
KA = k .* A;
G = bsxfun(@rdivide, KA, sum(KA, 2));
end